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Abstract 

A comparison is made between the Hagen-Poiseuille flow in rigid tubes and 
networks on one side and the time-independent one-dimensional Navier- 
Stokes flow in elastic tubes and networks on the other. Analytical rela- 
tions, a Poiseuille network flow model and two finite element Navier-Stokes 
one-dimensional flow models have been developed and used in this inves- 
tigation. The comparison highlights the differences between Poiseuille and 
one-dimensional Navier-Stokes flow models which may have been unjustifi- 
ably treated as equivalent in some studies. 

Keywords: fluid mechanics; Poiseuille; ID flow; Navier-Stokes; Newtonian 
fluid; rigid tube; distensible tube; network of tubes; porous media; hemody- 
namics. 


1 Introduction 

Poiseuille flow model has been widely used in the earth science and engineering 
studies for modeling and simulating the flow of Newtonian fluids in networks of 
rigid tubes which may represent a network of interconnected pipes for oil trans- 
portation or a simplified imitation of porous media. The model is naturally ex- 
tended to include Poiseuille-like flow of time- independent non-Newtonian fluids 
[1-7]. Although the ID Navier-Stokes flow model for elastic tubes and networks 
is the more popular [8-17] and normally is the more appropriate one for biologi- 
cal hemodynamic modeling, Poiseuille model has also been used in some studies 
for modeling and simulating blood flow in large vessels without accounting for the 
distensibility of the biological networks. 

One of the main differences between Poiseuille and ID flow models for single 
tube, which also affects the network flow since the individual tubes in the network 
are subject to the same flow principles as the stand-alone tubes, is that the flow 
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in the first model depends on the pressure difference while in the second model it 
depends on the actual pressure at the inlet and outlet [18] . This difference is mainly- 
based on the rigidity and distensibility of the flow ducts in these two types of tubes 
and networks. Another principal difference between Poiseuille and ID models for 
network flow, which reflects their complexity and practical relevance and originates 
from their single tube models, is that Poiseuille is linear and hence it is numerically 
solved in a single iteration, while the ID model is nonlinear and hence it requires an 
iterative process which may cause convergence instabilities leading to compromises 
associated with considerable numerical errors and approximations. 

Apart from the appropriateness of one of these models or the other for a given 
physical situation (a reason that dictates which model must be used in a specihc 
situation) there are certain practical advantages and disadvantages of Poiseuille 
and ID network flow models. In general Poiseuille is easier to implement, more- 
over it incurs a relatively low computational cost, normally of the order of N'^ of 
memory space where N is the number of network nodes, while the ID is more dif- 
hcult to implement with a high computational cost of the order of 4A^^ of memory 
space. This memory cost, associated with the previously-indicated solver iteration 
requirement, have obvious consequences on the speed of operation and overall per- 
formance. The high computational cost of the ID model in terms of memory space 
and CPU time can substantially increase by the demand of hne meshing and the 
use of higher order interpolation. 

On the other hand, the ID model gives more detailed picture as it depicts the 
flux and pressure helds over the spatial domain inside the tubes; opposite to the 
Poiseuille model which can only provide the average flow rate in the tubes and 
the pressure at the junctions. However, the pressure at the interior points of the 
tubes can be obtained directly for Poiseuille flow due to the linearity of the pressure 
held; a feature that cannot be replicated in the ID how due to the nonlinearity 
of the pressure held. Another advantage of the ID Navier-Stokes model is that it 
can simulate transient how as well as steady state how, while the Poiseuille model 
is basically time-independent and hence any time-dependent feature cannot be 
directly and dynamically replicated. However, consideration of transient ehects in 
Poiseuille model may be imitated indirectly by the generation of a sequence of time- 
independent how frames which represent snapshots of the overall time-dependent 
process. 
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2 Poiseuille Model for Single Tube and Network 


Poiseuille formula for the flow of laminar, incompressible, axi-symmetric, time- 
independent, fully-developed flow of Newtonian fluids in rigid cylindrical tubes 
assuming no-slip-at-wall [19] conditions is given by 




Tir^AP 

8fiL 


( 1 ) 


where Q is the volumetric flow rate, r is the tube radius, AP is the pressure drop 
along the tube, /i is the fluid dynamic viscosity and L is the tube length. This 
equation can be derived by several methods [1, 4, 5, 20]; most of which are based 
on the use of Navier-Stokes equation or one of its subsidiaries. 

There are several methods for modeling and implementing Poiseuille flow in a 
network of rigid tubes. The methods are mainly based on imposing a Poiseuille flow 
condition dehned by a constant conductance, as given by Equation 1, on each tube 
in the network associated with a conservation of flow condition on all junctions 
of the network, as well as boundary conditions on all the boundary nodes. In 
the following we briefly describe one of these methods which may be the most 
straightforward one and is based on imposing pressure boundary conditions. 

For a network with N nodes, which include all the internal junctions as well 
as the boundary nodes, the following system of simultaneous linear equations is 
formed 


CP = Q (2) 

where C is the tubes conductance matrix with dimensions of x iV, P is the 
pressure column vector with dimensions of x 1 and Q is primarily the total flow 
column vector with dimensions of A^ x 1. The rows of C basically represent the 
relations describing the net volumetric flow rate on the nodes in the network, which, 
apart from the boundary nodes, sums up to zero due to the incompressibility of 
fluid with the absence of sources and sinks inside, while the columns represent the 
nodes of the network. Similarly, the rows of the P vector represent the pressure at 
the N nodes while the rows of the Q vector primarily represent the net flow at the 
nodes. The pressure drop across each tube is then split into two parts corresponding 
to the two nodes involved at the two ends of that tube. All the entries in each row 
of the conductance matrix are set to zero except the ones related to the tubes 
that are connected to the node represented by that row. By a proper choice of 
flow signs, the total flow at each internal junction will be added to zero to satisfy 
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the continuity condition and hence the only sources and sinks for the flow are the 
boundary nodes. For the rows corresponding to the boundary nodes, the boundary 
pressure conditions are imposed by setting the diagonal entry of the conductance 
matrix to unity while all the other entries are set to zero with the corresponding 
entry of the Q vector being set to the given pressure boundary value of that node. 

Poiseuille model for the networks can be validated by testing the satisfaction of 
Poiseuille formula (Equation 1) on all the tubes of the network given the obtained 
pressure at the two end nodes of each tube, plus the satisfaction of the boundary 
conditions on the boundary nodes and the continuity of flow at each junction node 
which requires that the net flow in each internal node should be added to zero. 
A consequence of this is that the outflow (sum of volumetric flow rate at outlet 
boundary nodes) should be equal to the inflow (sum of volumetric flow rate at inlet 
boundary nodes). 


3 ID Model for Single Tube and Network 


The one-dimensional model for the time- independent flow of Newtonian fluids in 
elastic tubes is derived from the following Navier-Stokes system, which is based on 
the mass and momentum conservation principles 
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In these equations, A is the tube cross sectional area, Q is the volumetric flow 
rate, z is the axial coordinate along the tube, L is the length of tube, a (= 
with u and u being the fluid local and mean axial speed respectively) is the correc- 
tion factor for momentum flux, p is the fluid mass density, p is the local pressure, 
and K is the viscosity friction coefficient usually given hj k = 2nau/{a — 1) with 
u being the fluid kinematic viscosity dehned as the ratio of the dynamic viscosity 
p to the mass density. The model is based on the assumption of a laminar fully- 
developed flow of incompressible fluids with a tube having a pressure-dependent 
cross sectional area. 

The relation between the pressure and flow rate in elastic tubes is dependent on 
the constitutive relation that links the local pressure at a certain point along the 
tube axis to the cross sectional area of the tube at that point. In reference [18] the 
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following two flow relations have been derived based on two different pressure-area 
constitutive relations 


Q 
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and 
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In these equations, Ao is the tube unstretched cross sectional area at reference 
pressure, while Ain and Aou are the cross sectional area at the inlet and outlet re- 
spectively. Equations 5 and 6 are based on the following pressure-area constitutive 
relations respectively 


where, in the last two equations, p is the actual pressure relative to the reference 
pressure with which the unstretched area in defined, and 7 and (3 are the propor- 
tionality coefficients that control the stiffness of the tube. 

On multiplying the mass and momentum conservation equations by weight func- 
tions and integrating over the solution domain the weak form of the Navier-Stokes 
flow system can be obtained. This weak form, with suitable boundary conditions, 
can then be used as a basis for finite element implementation in conjunction with 
an iterative scheme such as Newton-Raphson method. A detailed account about 
the hnite element formulation is given in references [18, 21]. 

The system can also be extended to a network of elastic tubes by imposing 
suitable pressure or flux boundary conditions on all the boundary nodes, and com- 
patibility and matching conditions on all the internal junctions, where the latter 
conditions are derived from Riemann’s method of characteristics, and mass and 
energy conservation [16, 21, 22]. 

The time-independent hnite element solution of the ID model on a network 
is validated by fulhlling the boundary and matching conditions plus the time- 


p = 7 (A - Ao) 


(7) 


and 



( 8 ) 
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independent analytical solution, as given by Equations 5 and 6, on each individual 
tube. The validation of the time- independent hnite element solution of the ID 
model on a single tube is trivial as the numeric solution can be compared directly 
to the analytic solution with the inspection of the boundary conditions. 

4 Comparing Poiseuille and ID Models for Sin- 
gle Tube 

Thorough comparisons between rigid Poiseuille and elastic ID Navier-Stokes mod- 
els have been carried out for single tube as part of this study to investigate the 
effect of various parameters of these two models on the flow behavior. The results 
of a sample of these comparisons, which represent the two previously mentioned 
pressure-area constitutive relations, are plotted in Figures 1, 2, 3 and 4. The 
Poiseuille results shown in these hgures are obtained from the analytical solution 
(Equation 1) while the ID results are obtained from the analytical expressions 
(Equation 5 for the hrst p-A relation and Equation 6 for the second p-A relation) 
with an endorsement by a hnite element numeric implementation [18]. 

As seen, the investigated how, huid and tube parameters; which include a, /?, 7, 
/i, pressure regime, p, and tube size; have signihcant ehects on the how conduct of 
these two models and hence the results are highly dependent on these parameters. 
It is noteworthy that all these reported results represent qualitative demonstration 
and hence may not rehect a general trend due to the ehect of other parameters 
which are held constant to investigate the particular dependency. The variation 
of these parameters is expected to ahect the apparent trend in general. We also 
do not report the values of the other parameters for each one of these cases in 
detail due to the generality of the current study and the qualitative nature of this 
demonstration, as well as space limitation and avoiding unnecessary repetition. 
However, in all these simulations, typical and representative values have been used 
for the parameters related to the flow, fluid and flow paths, unless stated otherwise. 
This general conduct has also been followed in the forthcoming investigation of the 
flow in networks. 
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(e) fi = 0.05 Pa.s (f) ^ = 0.1 Pa.s 


Figure 1: The effect of a, 7 and /i on the difference between Poiseuille model 
for rigid tube and the ID Navier-Stokes model for elastic tube using the first p-A 
relation. The planar surface belongs to Poiseuille while the curved surface belongs 
to the ID. The horizontal axes represent the pressure at the two ends of the tube 
where the inlet is taken where the largest pressure occurs and hence the two surfaces 
are symmetric with respect to the diagonal zero-ffow-rate line of equal pressures. 
The pressures are in Pa while the flow rates are in m^/s. 
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(a) Pressure range: 0-1000 Pa 


(b) Pressure range: 0-2000 Pa 





(e) r = 0.1 m, L = 1.0 m 


(f) r = 0.01 m, L = 0.1 m 


Figure 2: The effect of pressure range, p and tube size on the difference between 
Poiseuille model for rigid tube and the ID Navier-Stokes model for elastic tube 
using the first p-A relation. All the other aspects are as in Figure 1. 
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(a) a = 1.1 


(b) a = 1.333 



(c) /3 = 10^ Pa.m 


(d) /3 = 10^ Pa.m 



1000 


(e) fi = 0.05 Pa.s 



(f) pL = 0.1 Pa.s 


Figure 3: The effect of a, /3 and /i on the difference between Poiseuille model for 
rigid tube and the ID Navier-Stokes model for elastic tube using the second p-A 
relation. All the other aspects are as in Figure 1. 
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(a) Pressure range: 0-1000 Pa 


(b) Pressure range: 0-2000 Pa 




1000 



1000 


(e) r = 0.1 m, L = 1.0 m 


(f) r = 0.05 m, L = 0.5 m 


Figure 4: The effect of pressure range, p and tube size on the difference between 
Poiseuille model for rigid tube and the ID Navier-Stokes model for elastic tube 
using the second p-A relation. All the other aspects are as in Figure 1. 
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5 Comparing Poiseuille and ID Models for Net- 
work 

Both p-A models have been implemented in a finite element network fiow computer 
code as part of this study. However, due to instabilities and convergence difficul- 
ties encountered in some cases of the first model which may affect the reliability 
and generality of the results, we will only present and analyze the results of the 
second model. This will not affect the generality of this study since a single p-A 
model is sufficient to highlight the main issues and draw the conclusions; moreover 
the generic trends observed from both models are similar in general. In fact the 
second model is the one which is widely used by the elastic network fiow modelers, 
possibly due to its numerical stability and reliability. The effect of the p-A model 
on the Poiseuille versus ID comparison may deserve a study by itself. It should be 
remarked that all the network results reported in the current paper, like the single 
tube results, have been subjected to rigorous checks using the validation criteria 
outlined in sections 2 and 3. 

With regard to the type of networks used in this investigation for Poiseuille and 
ID fiow simulations, we used fractal-type networks generated by a computer code. 
The main feature of these networks, which differ in the number of generations 
and consequently the number of elements which ranges between a few tubes to 
several thousands, is that they have a fixed branching angle with a constant length 
to radius ratio; moreover the branching radius transition from one generation of 
tubes to the next generation is subjected to a Murray-type rule [21], i.e. 

n 

T = E d (9) 

i 

where is the radius of the mother tube, is the radius of the ith. daughter tube, 
n is the number of daughter tubes, and 5 is a constant index. Various networks 
with different branching angle, length-to-radius ratio, Murray’s index, and number 
of generations have been used in our investigation. However, most of our networks 
were generated with n = 2 and 6 between 2 and 3 with equal-size daughter tubes. 
A sample of these networks with different number of generations are presented in 
Figure 5. The reason for using highly regular and symmetric fractal networks is that 
they, with their simplicity and regularity, reduce the effect of factors related to the 
complex structure. With the use of these simple fractal networks the fiow results 
will essentially reflect the correlation of the fiow with the varied parameters in a 
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simple manner. The nse of networks with complex morphology will only obscure 
the results and complicate the analysis due to the involvement of factors related to 
the complex geometry and topology of the network. 

Using these fractal networks in conjunction with the second p-A model for the 
ID Navier- Stokes, a number of Poiseuille and ID simulations have been conducted 
using typical network, fluid and flow parameters. In these simulations a typical 
case in which Poiseuille and ID produce virtually identical results with a = 1.333, 
= 236 Pa.m, p = 0.0035 Pa.s, Pi = 2000 Pa, Po = 0 Pa, and p = 1060 kg/m^ has 
been used as a reference case. To investigate the effect of the parameters of concern 
{a, f3, p, p, p, and size) a single parameter was varied at a time from this reference 
case and a comparison between Poiseuille and ID was made to identify the effect of 
that parameter. The method of comparison is based on computing the ratio of the 
volumetric flow rate between Poiseuille and ID flow for each tube and plotting this 
ratio versus tube indices, as seen in Figure 6. Although using the ratio of pressure 
at each junction node is sensible for making the comparison it has not been used 
here due to the fact that the pressure at each node for the ID model is tube- 
dependent because each tube at a specihc junction has its own pressure due to the 
involvement of Bernoulli energy conservation principle as a matching condition. As 
for the boundary conditions which were employed in these simulations, a single inlet 
node belonging to the largest single tube was used to impose an inlet Dirichlet-type 
pressure boundary condition while all the other boundaries at the other end, which 
belong to the smallest tubes representing the last generation of the network, were 
subjected to zero-pressure boundary conditions. However, there is one exceptional 
case where the outlet pressure was set to a non-zero value to investigate the effect 
of pressure limits, as seen in Figure 6(e). 

On inspecting Figure 6 it can be seen that while a and have a signihcant 
effect on the ID flow model as compared to Poiseuille, the other parameters have 
either moderate or negligible effect. However, this may not be true in general due 
to the limitation of this study and the number of cases investigated. The use of 
networks with more complex morphology is expected to introduce other sources of 
discrepancy between the two models and exacerbate the difference between them. 
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(c) 10-generation (d) 15-generation 


Figure 5: A sample of fractal networks used in the current investigation with 
different number of generations as well as other branching parameters. 
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Figure 6: Comparing the effect of various parameters on Poiseuille and ID Navier- 
Stokes network flow models using the second p-A constitutive relation. The |/-axis 
represents the ratio of Poiseuille flow rate, Qp, to ID flow rate, Qn- 
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6 Conclusions 


The purpose of this study is to compare Poiseuille model for rigid tubes and net- 
works with the time-independent ID Navier-Stokes model for elastic tubes and 
networks based on investigating the effect of the parameters of these two models 
related to the flow, fluid and flow ducts. The main conclusion of the current inves- 
tigation is that Poiseuille and ID models could produce very different results and 
hence they should not be used interchangeably as it may happened in some stud- 
ies. The use of one model or the other should be based on the merit of that model 
and its suitability to capture the physical reality and similar objective consider- 
ations, not on convenience and pragmatic factors. Moreover, the results of these 
models should be assessed relying on independent metrics such as consistency and 
compliance with experimental observations. 

Apart from the main theme of this investigation, it has been observed that 
for the investigated cases of single tube, all the investigated parameters; which 
include a, /3, 7 , p, p, pressure range, and tube size; have signihcant effects. With 
regard to the networks, limited in this study to those with fractal character, it has 
been observed that a and f3 have the most sizeable impact. However, the impact 
of each one of the investigated parameters can be affected by the involvement of 
other factors such as network topology and geometry. Moreover, other factors are 
expected to play a signihcant role in networks with more complex morphology. 

Other effects related, for example, to the converging-diverging nature of the 
how ducts and non-Newtonian rheology [23-25], are not considered in this study. 
These and all other factors should contribute to the complexity of the situation and 
the departure of the results of Poiseuille and Poiseuille-like models for rigid tubes 
and networks from the results of the ID models for distensible tubes and networks. 

The comparison presented in the current paper is very general and lacks thor- 
oughness due to specihcity of purpose and space limitation. The ehect of each 
one of the investigated parameters, as well as many other aspects not touched in 
this study such as the ehect of morphology of the networks and their statistical 
distribution, deserve a study on its own. 


15 


Nomenclature 


a correction factor for axial momentum flux 
(3 stiffness parameter in the second pressure-area relation 
7 stiffness parameter in the hrst pressure-area relation 
5 Murray’s law index 

K viscosity friction coefficient 

jj, fluid dynamic viscosity 

V fluid kinematic viscosity 

p fluid mass density 

A tube cross sectional area 

Ain tube cross sectional area at inlet 

Ao tube cross sectional area at reference pressure 

Aou tube cross sectional area at outlet 

C conductance matrix 

L tube length 

n number of daughter tubes in Murray’s law 

p local pressure 

P pressure column vector 

Pi inlet pressure 

Po outlet pressure 

AP pressure drop along the tube 

Q volumetric flow rate 

Q total flow column vector 

Qn volumetric flow rate of ID Navier-Stokes model 
Q p volumetric flow rate of Poiseuille model 
r tube radius 

Td radius of daughter tube 

Tm radius of mother tube 

u local axial speed of fluid 

u mean axial speed of fluid 

tube axial coordinate 
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